Exploration and identification of anoikis-related genes in polycythemia vera

Background: Polycythemia Vera (PV) is a type of typical Myeloproliferative Neoplasms (MPNs) characterized with excessive erythropoiesis and thrombosis. Anoikis is a special programmed cell death mode induced by the adhesion disorder between cells and extracellular matrix (ECM) or adjacent cells facilitating cancer metastasis. However, few studies have focused on the role of anoikis in PV, especially on the development of PV. Methods: The microarray and RNA-seq results were screened from the Gene Expression Omnibus (GEO) database and the anoikis-related genes (ARGs) were downloaded from Genecards. The functional enrichment analysis of intersecting differentially expressed genes (DEGs) and protein-protein interaction (PPI) network analysis were performed to discover hub genes. The hub genes expression was tested in the training (GSE136335) and validation cohort (GSE145802), and RT-qPCR was performed to verify the gene expression in PV mice. Results: In the training GSE136335, a total of 1,195 DEGs was obtained from Myeloproliferative Neoplasm (MPN) patients compared with controls, among which 58 were anoikis-related DEGs. The significant enrichment of the apoptosis and cell adhesion pathways (i.e., cadherin binding) were shown in functional enrichment analysis. The PPI network was conducted to identify top five hub genes (CASP3, CYCS, HIF1A, IL1B, MCL1). The expression of CASP3 and IL1B were significantly upregulated both in validation cohort and PV mice and downregulated after treatment, suggesting that CASP3 and IL1B could be important indicators for disease surveillance. Conclusion: Our research revealed a relationship between anoikis and PV for the first time by combined analysis of gene level, protein interaction and functional enrichment, allowing novel insights into mechanisms of PV. Moreover, CASP3 and IL1B may become promising indicators of PV development and treatment.

Anoikis, a type of programmed apoptosis, is induced by incorrect attachment between cells and extracellular matrix (ECM) or adjacent cells, playing important role in eradicating misplaced or dislodged cells to maintain tissue homeostasis (Han et al., 2021). Some studies suggested anoikis occurred through intrinsic/extrinsic pathways and could be triggered by several intracellular signals such as DNA damage and endoplasmic reticulum stress (Li et al., 2019;Zhi et al., 2022). Anoikis is essential for cancer cell survival through the detachment from ECM, and anoikis resistance is a key factor in cancer invasion and metastasis (Berezovskaya et al., 2005;Kim et al., 2012;Kakavandi et al., 2018;Xiao et al., 2019;Yu et al., 2022;Zhou et al., 2022).
Clonal cell proliferation leads to increased risk of thrombosis and threatens the life of PV and ET patients (Barbui et al., 2013). Whereas, complications of thrombosis in myelofibrosis (MF) are as frequent as in ET but less than in PV (Kc et al., 2017). Recently, the abnormal expression of adhesion molecules like Lu/BCAM and CD147 was reported in PV patients, which was reduced under interferon-α treatment but enhanced by hydroxyurea (Brusson et al., 2018). These results indicated the important role of cell adhesion in PV and PV treatment. However, few studies had systematically evaluated the effect of anoikis in hematological cancers especially MPNs. Therefore, in this study, we focused on exploring the expression levels of anoikis-related genes in the occurrence, progression and treatment of MPNs, as well as exploring the key anoikis-related genes by means of bioinformatics analysis.

Data collection
Key word "myeloproliferative neoplasms" was used to search gene expression profiles of MPNs in the Gene Expression Omnibus (GEO, http://www.ncbi.nlm.nih.gov/geo/) data portal. The following criteria were used to filter the obtained database: 1. The gene expression profiling must cover all types of chronic and progressed MPNs (PV, ET, PMF) and controls. 2. The cell used for sequencing should be CD34 + MNCs (mononuclear cell) from peripheral blood (PB) or bone marrow (BM). 3. The processed data or raw data from datasets must be provided thus could be available to reanalyze. 4. The GEO dataset numbered GSE136335, GSE145802 were selected. Further, anoikis-related genes (ARGs) were downloaded from the GeneCard database (https://www. genecards.org/).

Design of this study
The study design was shown in Figure 1. Firstly, differentially expressed genes (DEGs) were screened in MPNs and MPN sub types including PV, ET, PMF and post-PV myelofibrosis (pPVMF), as well as the intersection of the DEGs and ARGs was analyzed. And then functional enrichment GO and KEGG analyses were performed. Lastly, PPI network analysis was conducted and the top 5 hub genes were found in the network. In addition, the corresponding hub genes in different validation cohorts were validated and the hub genes of PV in our mice model were detected.

Data processing and screening of anoikisrelated DEGs
Differentially expressed genes (DEGs) between patients with MPNs and healthy controls were identified in GSE136335 dataset via the "limma" R package (Gentleman, 2005), the cutoff for adjustment: p-value <0.05 and FC (fold change) > 1.5. The "ggplot2" R package (Ito and Murphy, 2013) was used to visualize the DEGs into a volcano map, while the "pheatmap" R package (Szekely and Rizzo, 2005) was used to cluster the anoikisrelated DEGs. The intersection between DEGs and ARGs was visualized by the Veen plot.
Functional enrichment analysis and protein-protein interaction network analysis of anoikis-related DEGs.
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of anoikis-related DEGs were performed by using "clusterProfiler" R package . GO analysis is important to explore the biological functions, consisting of three subtypes, including biological process (BP), cellular component (CC), and molecular function (MF) (Ashburner et al., 2000). KEGG analysis is used to explore potential pathways (Kanehisa and Goto, 2000) The STRING online website (https://string-db.org/) was to conduct the protein-protein interaction (PPI) network analysis of anoikisrelated DEGs. The results of STRING were imported into Cytoscape (version 3.7.2), and the hub genes was extracted by using the cytoHubba plugin. Finally, top five genes were selected with the highest scores (Degree algorithm) in the network as hub genes.

Mouse model
Animal experiments were approved by the Central South University, Changsha, Institutional Animal Care and Use committee. Animals were maintained a 12:12 h light/dark cycle Frontiers in Genetics frontiersin.org with a temperature-controlled specific pathogen free (SPF) environment. 8-wk-old male mice were used in this study. For analyses in the Jak2V617F bone marrow transplantation (BMT) model, experiments were performed in strict adherence to China laws for animal welfare. This study was approved by the Ethics Committee of The Second Xiangya Hospital of Central South University (20220247). For gene expression analysis in the conditional Jak2V617F knock-in model of PV, bone marrow (BM) from primary Vavcre Jak2V617F (Vavcre + Jak2 VF/+ ) mice (Tiedt et al., 2008;Mullally et al., 2010) (Jackson Lab (United States) was transplanted into lethally irradiated C57BL/6 mice (Slac Laboratory Animal Inc. Shanghai, China) (PV group), and the recipients progressed to PV before being used in the MEPs (megakaryocyte and erythroid progenitors) flow sorting experiment. As competitor transplantation, bone marrow from C57BL/6 recipient mice was transplanted into lethally irradiated C57BL/6 recipients (WT group).
6-8 weeks after BMT, BM cells of PV and WT recipients were harvested from bilateral femurs and tibias by crushing bones with mortar and pestle using PBS (Servicebio). Cells were filtered through 70 μm nylon mesh to obtain a single-cell suspension. Red blood cells were depleted by using erythrocyte lysis buffer (ACK buffer, Invitrogen).
For FACS analysis and cell sorting: EasySep ™ Mouse Hematopoietic Progenitor Cell Isolation Kit (Catalog #19856, STEMCELL) was used to enrich the Lin − cells with biotinylated antibodies recognizing specific cell surface markers. The cell surface was stained using antibodies as shown in Supplementary Table S2: (Mead et al., 2013;Rao et al., 2019). 80,000-120,000 MEPs were isolated from PV and WT recipient mice by FACS sorting and used to further gene expression analysis. Cells were analyzed on a Fortessa Flow Cytometer and sorted on a FACSAria-II cell sorter (BD biosciences). Data was analyzed either using FACS Diva software or FlowJo (version 10.8.1) software (Treestar).

RNA extraction and quantitative realtime PCR
After extracting total RNA of flow-sorted MEP cells using TRIzol (Invitrogen, United States) according to the manufacturer's protocol, reverse transcription high-capacity cDNA reverse transcription kit (AG, AG11728, Changsha, China) was performed. Lastly qRT-PCR using SYB Green Master quantitative PCR system (YEASEN, Shanghai, China) was conducted. Relative expression levels was calculated by using the 2 -△△CT method. The sequences of the primers were listed in Supplementary Table S3.

GEO information
The study flow chart was shown in Figure 1. The two GEO datasets numbered GSE136335 and GSE145802 were selected to study according to the criteria previously set. Information and characteristics of the three datasets were summarized in Table1. The GSE136335 was confirmed as a discovery cohort with GSE145802 paired as a validated cohort.

Identification of anoikis-related DEGs between MPNs and controls
1195 DEGs were identified by screening GSE136335, with 557 downregulated genes (Figure 2A, blue dot) and 638 upregulated genes (Figure 2A, red dot). 795 anoikis-related genes were downloaded from the "GeneCard" database (Supplementary Table S1). Veen diagram ( Figure 2B); the left region (1,137 genes) refered to unique DEGs in GSE136335; the middle intersection region (58 genes) refers to anoikis-related DEGs; and the right region (736 genes) refered to unique anoikis genes. The heat map in Figure 2C showed the standardized expression of anoikis-related DEGs (36 upregulated and 22 downregulated). The PPI network ( Figure 2D) containing 54 nodes and 144 edges was obtained. Figure 3 showed the relative expression profiles of 58 anoikis-related DEGs in the MPNs ( Figure 3A) and MPNsubtypes ( Figure 3B).
In GO-BP analysis ( Figure 4A), the major pathways were also apoptosis or oxidative stress pathways. In GO-CC analysis ( Figure 4B), collagen-containing extracellular matrix was significantly enriched. The results of enrichment analysis in GO-MF ( Figure 4C) were mainly the activation of cell adhesion-related pathway like cadherin binding. In the KEGG analysis ( Figure 4D), the top 10 enriched pathways were mainly apoptosis, FoxO signaling pathway, TNF signaling pathway, etc.

Identification of anoikis-related DEGs between PV and controls
925 DEGs were identified in PV patients compared with controls, with 513 downregulated genes and 412 upregulated genes ( Figure 5A). Intersection analysis showed 36 anoikisrelated DEGs ( Figure 5B). The standardized expression of anoikis-related DEGs (17 upregulated and 19 downregulated) was shown in the heat map ( Figure 5C). The results of enrichment analysis in GO and KEGG ( Figure 5D; Supplementary Figure S1A) were also mainly focused on the apoptosis and activation of cell adhesion-related pathway like integrin binding. PPI network ( Figure 6A) was constructed to select top 5 hub genes (CASP3, IL1B, HIF1A, CYCS, MCL1) ( Figure 6B) according to the Degree algorithm. The expression levels of hub genes in training set were shown ( Figure 6C).

FIGURE 3 Expression of 58 anoikis-related DEGs in MPN (A) and subtype clusters (B).
Frontiers in Genetics frontiersin.org Validation analysis was conducted by using another dataset GSE145802. The expression levels of anoikis related DEGs in PV, CASP3, IL1B, HIF1A were significantly upregulated and CYCS was obviously downregulated, while the expression of MCL1 was downward trend in PV patients compared with normal controls ( Figure 7A). However, the expression levels of anoikis related DEGs in PV-treated, only HIF1A was upregulate ( Figure 7B).

Expression difference of hub genes in PV mice model
The expression of five PV hub genes was verified with our experimental data. PV mice model with Jak2V617F mutation (PV) and the corresponding WT mice (WT) were constructed by bone marrow transplantation (BMT) assay ( Figure 8A). The spleen weight ( Figure 8B) and hematocrit (HCT) ( Figure 8C) of the recipient mice were measured 6-8 weeks after BMT, and the results indicated PV models successfully were constructed. PV and WT MEP cells were obtained by Mouse Hematopoietic Progenitor Cell Isolation Kit combined with FACS sorting (List of flow cytometry antibodies can be found in Supplementary Table  S2) following the strict experimental procedures ( Figure 8D). The relative mRNA expression (The sequences of the primers are listed in Supplementary Table S3) of Casp3 and Il1b were upregulated in MEPs of PV mice compared with WT mice ( Figure 8E). However, other hub genes had no difference between these two group mice.

Discussion
Myeloproliferative neoplasms (MPNs) are myeloid malignancies characterized by clonal expansion of one or more differentiated myeloid lineages (Dameshek, 1951). Anoikis, a particular form of programmed cell death, was reported to perform important biological function in preventing attachment Frontiers in Genetics frontiersin.org or grow of dysplastic cell for tissue homeostasis and development (Bakir et al., 2020;Li et al., 2021). Anoikis plays a crucial role in tumor invasion, migration, distant metastasis and drug resistance (Jin et al., 2017;Jiang et al., 2020;Kim et al., 2021). In our study, a new view of differential expression anoikisrelated gene sets in MPN patients and controls was provided. MPNs are the malignancies caused by dysfunctional HSC, thus CD34 + PB/BM MNCs was chosen for bioinformatics analysis. The intersection of MPN-DEGs and ARGs revealed 58 anoikisrelated genes were specifically associated with MPNs. Gene set enrichment analysis showed cell adhesion and apoptosis-related pathways were fully enriched to activate the anoikis pathway. However, in MPNs, circulating cells with disease driver mutations (such as JAK2V617F), including red cells, leukocytes, platelets, as well as some vascular endothelial cells, were abnormal (Hasselbalch et al., 2021). When activated in MPNs, these abnormal cells interacted with each other and created proadhesive and prothrombotic environment in the circulation to predispose patients to thrombosis and occlusive complications (Hasselbalch et al., 2021).
PV is characterized by mutations in JAK2 exon 14 or 12 in hematopoietic stem cells, leading to erythrocytosis and systematic symptoms (Greenfield et al., 2021). PV mainly reduces the quality of life and shortens the life span of patients due to the risk of vascular events and transformation to myelofibrosis, myelodysplastic syndrome, and acute myeloid leukemia (Tefferi et al., 2013;FIGURE 5 Landscape of Anoikis-related Genes in Polycythemia Vera. (A)Volcano plot of genes differentially expressed between Polycythemia Vera (PV) patients and healthy controls in the GSE136335 dataset. Blue nodes represent downregulation in PV; red nodes represent upregulation; and grey nodes represent no significant difference from controls. (B) Intersection of differentially expressed genes (DEGs) in the GSE136335 and anoikis genes. The count on the left (889 genes) refers to DEGs unique to GSE136335; the count in the middle (36 genes) refers to anoikis-related DEGs; and the count on the right (758 genes) refers to unique to anoikis genes. (C) Heat map of 36 anoikis-related DEGs. (D) Top 10 GO biological process pathway.

Frontiers in Genetics
frontiersin.org

FIGURE 6
The protein-protein interaction network of anoikis-related DEGs in Polycythemia Vera. (A) PPI network for 36 anoikis-related DEGs; (B) Subnetwork of hub genes from the PPI network. Node color reflects the degree of connectivity (red color represents a higher degree, and yellow color represents a lower degree). (C) The expression difference of five hub genes in the training set (GSE136335). Frontiers in Genetics frontiersin.org 08 Tiribelli et al., 2015;Frederiksen et al., 2019). Our results found that 36 DEGs in PV were associated with anoikis. Combined with the results of gene expression and gene enrichment, CD34 + mononuclear cells in PV deflected stronger adhesion effect than apoptosis, in accordance with the characteristics of PV prone to thromboembolism. Five hub genes (CASP3, IL1B, HIF1A, CYCS, MCL1) was selected with highest score through PPI network combined further analysis. The pro-apoptotic gene CASP3 was expressed higher in PV-CD34 + cells in the training and validation cohorts and in MEP cells in our PV mouse model, contrary to previous report (Steidl et al., 2005). The large variation may root in PV patients of different sources. The level of proinflammatory factor IL1B was upregulated in PV-CD34 + cells in training and validation cohorts and in MEP cells in our PV mouse model, consistent with previous report (Vaidya et al., 2012). Baumeister, Julian et al. reported higher protein level of HIF1α in 32D Jak2V617F-expressing cells than Jak2WT, and identified HIF-1 target genes involved in the Warburg effect as a possible underlying mechanism in JAK2V617F-positive MPN (Baumeister et al., 2020). The mRNA level of HIF1A was higher in PV-CD34 + cells in validation cohort. MCL1, a typical anti-apoptotic gene, plays a crucial role in the survival of HSCs and early progenitors (Opferman et al., 2005). Importantly, depletion of Mcl-1 was able to disrupt the viability of JAK2V617F mutant cells and sensitized mutant cells to JAK2 inhibitor therapy (Rubert et al., 2011). CYCS, as an apoptotic protein, is located between the inner and outer mitochondrial membranes (Ow et al., 2008;Shakeri et al., 2017) and is released into the cytoplasm in response to cell death signals and activated by pro-apoptotic Bcl2 family proteins (Wolter et al., 1997;Kennedy et al., 1999). The expressions of apoptotic-related genes MCL1 and CYCS decreased in the PV-CD34 + cells of the training and verification cohorts, while no difference was found in the PV mice MEP cells. CASP3 and IL1B were significantly increased in PV and decreased after treatment in validation cohort and expected to become indicators in PV development and treatment.
The inconformity between our results and previously reported may arise from following aspects: Firstly, the analyzed data was from peripheral blood or bone marrow mononuclear cells of MPN patients, excluding erythrocytes and granulocytes. Secondly, some vascular endothelial cells with driver mutations were also excluded. Thirdly, MPN in our study was a group of diseases rather than a single disease, and the degree of cell adhesion and apoptosis varies among subtypes. Finally, our results suggested that CD34 + mononuclear cells in bone marrow or peripheral blood contributed significantly less to thromboembolism than other types of cells.
Although MPNs are a myeloid proliferative disease, our research firstly revealed the critical role of anoikis-related genes in MPNs. CASP3, CYCS, HIF1A, IL1B, MCL1 are the anoikis-related hub genes in PV and CASP3 and IL1B may become important indicators of PV development and treatment. And our research may provide new insight into mechanisms of PV.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

Ethics statement
The animal study was reviewed and approved by the Ethics Committee of The Second Xiangya Hospital of Central South University.

Author contributions
WA, LX, and TD designed the study. WA and LX collected the data. WA and WH completed the experimental validation, including bone marrow transplantation, flow sorting and RNA extraction and quantitative realtime PCR. LX analyzed the data and WA drafted the manuscript. LX, TD, WH, and YT made critical revisions to the manuscript for important intellectual content. TD, HP, and GZ provided guidance on statistical analysis and presentation of results. All authors approved the final manuscript and agreed to submit for publication.

Funding
The National Natural Science Foundation of China, Key Program (82130024), the Science and Technology Innovation Program of Hunan Province (2020RC4009) and the Project of Innovation-Driven Plan of Central South University (2020CX015).